A cascade of magnetic field induced spin transitions in LaCo0 3 
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We present magnetization and magnetostriction studies of the insulating perovskite LaCo03 in 
magnetic fields approaching 100 T. In marked contrast with expectations from single-ion models, 
the data reveal two distinct first-order spin transitions and well-defined magnetization plateaux. 
The magnetization at the higher plateau is only about half the saturation value expected for spin- 
1 Co 3+ ions. These findings strongly suggest collective behavior induced by strong interactions 
between different electronic - and therefore spin - configurations of Co 3+ ions. We propose a model 
of these interactions that predicts crystalline spin textures and a cascade of four magnetic phase 
transitions at high fields, of which the first two account for the experimental data. 



PACS numbers: 

Spin state crossovers induced by changes in the elec- 
tronic configuration of transition metal ions can dramat- 
ically alter several materials properties [TJ [2] . One ex- 
ample is the pressure induced spin transition of ferric 
ions in magnesium silicate perovskite that occurs in the 
Earth's lower mantle [3j[4]. Apart from being a major 
constituent of the Earth's mantle, perovskites exhibit a 
variety of behaviors which include colossal magnetoresis- 
tance [5j [6] and high temperature superconductivity [7] . 
For certain insulating perovskites containing Fe or Co 
ions, the two lowest energy multiplets of the 3d electrons 
can be separated by a small energy gap, A, owing to 
competition between Hund's coupling and crystal field 
splitting. A external magnetic field, H, can change the 
electronic configuration if the lowest energy multiplet has 
lower spin than the first excited multiplet. In this situa- 
tion, structural transitions with unusually large magneto- 
elastic responses can be induced for H ~ A/gfis (where 
fiB is Bohr's magneton and g is the g-factor). 

The energy gap A separating low- and high-spin elec- 
tronic configurations of 3d ions in perovskites is typi- 
cally of order 10-1000 meV. Controlled switching between 
these electronic configurations with magnetic fields there- 
fore requires extremely large fields of order 100-10000 T 
(assuming g « 2). Therefore, such studies have proven 
difficult to date. The recent development of nondestruc- 
tive pulsed magnets with peak fields of 100 T and with 
long (ms-timescale) pulse durations now provides an op- 
portunity to realize field-tuned spin state transitions in 
compounds with A ~ 12 meV. Such studies are highly 
desirable, as they would provide an essential complement 
to pressure-induced spin cross-over phenomena that have 
been studied in a number of materials [8UTU]. Moreover, 
high-field studies allow one to test whether single-ion 
models suffice to describe spin-crossover phenomena in 



perovskites, or whether interactions between 3d ions are 
important. 

To this end we study the insulating perovskite cobaltite 
(LaCoOa, or LCO), whose octahedrally-coordinated 
Co 3+ (3d 6 ) ions are natural candidates to explore field- 
induced transitions of the electronic configuration and 
spin [TTJE]. A small energy gap A « 12 meV separat- 
ing the spin singlet ground state (6 electrons in the ti g 
orbitals) from the lowest energy magnetic configuration 
results from the competition between Hund's coupling 
and crystal field splitting [13] . While the Co 3+ ions are in 
their 5 = state at low temperatures, thermal activation 
to a magnetic 5 ^ state occurs above ~30 K, giving 
rise to a paramagnetic response. Although considerable 
work has been carried out on this thermally-induced spin 
crossover [T41JT6] . the spin value 5 of the lowest energy 
excited multiplet is still controversial. A single-ion model 
allows to explain the temperature dependence of the mag- 
netic susceptibility and the measured g-factor (g = 3.4) 
with an 5 — 1 triplet [15]. However, field-dependent 
magnetization studies have recently suggested that inter- 
actions between Co ions may also play a crucial role in 
LCO: a field-induced gap-closing study shows a jump in 
the magnetization M to a plateau value of ~ 0.5/i£ just 
above 60 T [17]. This value is much lower than the 2\±b 
(4/ig) expected from saturated 5 = 1 (5 = 2) Co 3+ ions 
within a single-ion model. From now on we will assume 
that 5 = 1, because this assumption leads to a simpler 
explanation of the observed plateaux. 

Very recently Platonov et al. [18 have reported the 
detection of several magnetic field-induced phase transi- 
tions in LCO single crystals at 4.2K. The samples are 
subjected to very high magnetic fields (« 500 T) pro- 
duced by rapid (« 15 /is) explosive compression of mag- 
netic flux. The results show that the magnetization 



per Co 3 + ion starts to rise significantly for H > 50T, 
where M/nco — 0.4/iB, reaching a plateau value of 
M/nco — 1.4/is at 140 T, and a maximum value of 
3.bfiB at 500 T. It is suggested that the observed plateau 
and other features may be linked to antiferromagnetic 
interactions. [T8] Smoothing of the first order transition 
near 60 T [17] is attributed to the rapid variation of 
the magnetic field. No firm conclusions on the Co 3+ 
spin state are presented, although it is suggested that 
the S = 2 state may be important at the highest fields 
reached. 

Here we measure the magnetization (M) of LCO in 
longer-duration non-destructive pulsed magnetic fields 
approaching 100 T [19]. In contrast to expectations 
from single-ion models, we observe multiple magnetiza- 
tion steps and plateaux, which indicate the relevance of 
inter-ion coupling. Moreover, magnetostriction studies 
complement the magnetization results, and reveal large 
lattice changes induced by a combination of two factors: 
a) the S=0 and S=l electronic configurations of the Co 3+ 
ions have different volume and b) the S=l configuration 
is Jahn- Teller active. Susceptibility measurements made 
at higher fields in a single turn pulsed field experiment re- 
veal that the second plateau persists at least up to 140 T. 
We propose a model that allows for collective behavior 
of the Co 3+ ions and reveals a rich phase diagram in 
which the plateaux originate from different stable crys- 
talline textures of S = 1 Co 3+ ions in a background of 
S = configurations. The orbital degree of freedom of 
the S = 1 Co ions also leads to orbital ordering that 
should coexist with the magnetic ordering. 

The single crystals of LCO are grown in a floating zone 
furnace at Argonne National Laboratory, with pieces be- 
ing cut and polished into either small tapered cylinders 
for the magnetization M measurements or 1 x 1 x 3mm 
bar shaped samples for the magnetostriction experi- 
ments. Magnetization curves are obtained as a func- 
tion of fioH, using an extraction magnetometer [20] . with 
pulsed fields up to 97 T at temperatures in the range 1.6- 
55 K. High precision magnetostriction measurements are 
made using a fiber-optic strain gauge [21] [22]. In all 
experiments, singe crystals are oriented with their [010] 
pseudo-cubic axes parallel to H. 

Figure [TJa) shows M versus hqH at T = 1.7 K re- 
vealing two transitions, each accompanied by increases 
AM ~ 0.5/iB/Co, and two plateau regions consistent 
with two distinct magnetic phases. The second plateau 
occurs at M ~ 0.9/i#/Co which is well below the mag- 
netization of 2 /i^/Co expected from saturated Co ions. 
Both transitions exhibit hysteresis in this pulsed field ex- 
periment. The location of the transition between rising 
and falling fields, as indicated by the blue and red arrows, 
points to a diffusionless, displacive transition type [23]- 
[26] . The inset in Fig. [l^a) for fields in the range 55- 
63 T, and with offsets for clarity, presents an expanded 
view of the first transition at temperatures between 1.6 
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FIG. 1: (Color online) (a) Magnetization of LaCoOa in units 
/jb/Co as a function of applied field at T = 1.7 K. The blue 
and red arrows indicate transitions for increasing and decreas- 
ing fields respectively. Two plateau regions can be seen. The 
inset shows an expanded view of the first transition between 
55 and 63 T for T = 1.6 (bottom), 4.0, 10.0, 20.0, 25.0, 27.0, 
28.0, 28.5, 29.0, 30.0, 40 and 55 (top) K. (b) Magnetostric- 
tion of LCO, measured using a high sensitivity optical strain 
gauge, as a function of applied field. The inset shows the 
susceptibility \ — dM/dH measured in a pulsed single turn 
coil (140 T in ~ 2 /xs). The two transitions are merged into 
a single peak in the very rapidly changing field. No further 
transitions are found below 140 T. 



K and 55 K. The hysteresis diminishes as T is raised. 
Once initiated, the transitions are observed to proceed 
rapidly, within ~ 5 /is, over the entire sample (length ~ 
4 mm). Figure [ljb) shows the magnetostriction mea- 
sured parallel to the field. The transition fields and 
AL/L plateaux closely match the magnetization behav- 
ior in Fig. [TJa). The red and blue arrows again represent 
rising and falling field transitions. The magnitude of the 
field-induced strain AL/L ~ 10 -3 exceeds that found in 
structurally related perovskites [27]. For the rapidly ris- 
ing fields, some rounding of the transitions is attributed 
to sample heating associated with latent heat released 
at the transitions, which occurs despite immersion in su- 
perfluid 4 He. The magnetostriction results, which show 
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FIG. 2: (Color online) Phase diagram for LCO based on the 
results shown in Fig. ma). Blue and red symbols denote rising 
and falling fields respectively, while black symbols connected 
by the solid curve represent the thermodynamic phase bound- 
ary obtained as an average of the rising and falling field values 
at a given T. SSC1 and SSC2 denote spin state crystalline 
one and two respectively. The inset shows low lying energy 
states as a function of applied field for Co 3+ in LCO based 
on a single ion description [TT] . 



sample expansion (dV '/ '8H)t,p > , are consistent with 
the Maxwell relations since (dM/dP) T ,H < for T < 100 
K [22] . The inset in Fig. [ljb) shows the magnetic suscep- 
tibility at 4 K measured in a pulsed single turn magnet 
which reaches ~ 140 T. 

Figure [2] shows a phase diagram for LCO based on the 
inset curves in Fig. [I] Besides the non-magnetic (NM) 
phase at low fields, there are two ordered magnetic phases 
denoted spin state crystalline 1 (SSC1) and 2 (SSC2) 
respectively. The inset depicts low-lying single-ion states 
for Co 3+ in LCO [15] HI]. 

The observed field induced phase transitions must in- 
volve interactions between the non-magnetic ions and the 
magnetic S z = +1 ionic configurations that are stabilized 
by the field. We note that the magnetic ions are both 
larger in volume as well as Jahn-Teller active [15] ITT] . 
Indeed, the S = +1 Co 3+ ion with a t\ g ^ g electron config- 
uration has an orbital degree of freedom since there is one 
electron shared between two e g orbitals. We assume that 
the energy gap (A) between the S = 1 and nonmagnetic 
configurations is large in comparison to the exchange in- 
teractions between different magnetic configurations of 
Co 3+ ions. We also assume that the gap to the S = 2 
magnetic configuration is much larger than 2A and con- 
sequently this configuration can be eliminated from the 
low-energy Hilbert space for fields up to 100 T. This as- 
sumption is supported by the consistency of the model 
predictions with experiment. Applying a magnetic field 
H lowers the energy of the doubly degenerate S z = +1 
states (the e g orbitals being d 3z 2_ r 2 and d x i_ y i) such 
that their energy becomes comparable to the S = non- 
magnetic state. A low-energy effective model is obtained 
by projecting the original Hamiltonian into the subspace 
of these three lowest-energy states. 



The orbital state of the S z = +1 ion is parameterized 
by an angle 0: \0) = cos((9/2)|3z 2 -r 2 )+sin((9/2)|x 2 -7/ 2 ), 
corresponding to a biaxial deformation of the local 06 
octahedra. The doublet orbital degeneracy manifested 
as a continuous symmetry for the angle is retained by 
linear Jahn-Teller (JT) interaction. In real compounds, 
however, the continuous degeneracy is lifted by the lat- 
tice anharmonicity and higher-order interactions [28]. In 
fact, only local elongations of the 06 octahedra are ob- 



served in almost all JT ions with e n electrons 



De- 



pending on the local elongation axis, the allowed orbitals 
are d 3x 2_ r 2, d 3y 2_ r 2 and d 3z 2_ r 2, corresponding to = 
and ±27r/3. The large nonlinear JT distortion also lifts 
the degeneracy of the t<i g levels and quenches the orbital 
degrees of freedom associated with the hole left in the t<i g 
manifold of the S z = 1 state. 

A single 5 = 1 ion embedded in the host matrix ac- 
quires additional energy via JT coupling with the lo- 
cal octahedron and exchange interaction with its non- 
magnetic neighbors, amounting to a slight modification 
(A) of the energy gap: 



^i = E( A -wff) 



(i) 



Here m = S£ = +1 such that rti = 1 for a magnetic con- 
figuration on the ion z, A is the renormalized spin state 
gap and gjisH is the Zeeman interaction. The single-ion 
physics described by Eq. ([T]) implies a single field-induced 
crossover when gjisH > A. To account for the observed 
multiple transitions, we need to include the interactions 
between the S z = +1 sites: 



n 2 



(2) 



The leading order coupling (Vij) is isotropic and repul- 
sive because it results from the increased volume of the 
S z = 1 ions relative to the nonmagnetic ones. The second 
interaction terms are orbital-dependent and determine 
the relative orbital orientations of the magnetic ions. [21] 
The three-state Potts variable Si a indicates the orbital 
states (d 3 i2_ r 2 with I = x : y, or z) of the S z = 1 ion. 

The origin of the isotropic interaction Vij can be 
understood by using the so called "sphere-in-the-hole" 
model [28]. The larger ionic size of a S = 1 ion in the 
host matrix of S = ions acts as an elastic impurity 
and creates a strain field which decays as ~ 1 /r 3 at large 
distances. A second impurity interacts with this strain 
field, giving rise to the elastic term in Equation Q. 

We find that the key aspects of the data — namely 
multiple phase transitions, incremental steps in the mag- 
netization, lattice expansions and metastability — can 
be captured by considering field-induced S z = +1 spin 
states in which we neglect orbital orientation-dependent 
terms. The Hamiltonian can then be mapped onto an 
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FIG. 3: (Color online) (a) The red curve gives the predicted 
normalized magnetization m versus /JoH behavior for LCO 
based on the interacting spin model. A cascade of four tran- 
sitions is shown at fields H c i to H C 4 as described in the text. 
The scaled green curve is based on the experimental results 
shown in Fig. [TJa). The right hand axis gives values of the 
pseudospin variable (a) used in the theoretical model, (b) 
Predicted spin state crystalline structures for each of the 
phases with fractional magnetization plateaux in Panel (a). 
The blue and (larger) red balls denote Co 3+ ions in the S = 
and S z = +1 states respectively. The intermediate spin Co 
ions form body-centered (BCC) and a face-centered (FCC) 
cubic structures in the m = 1 /4 and m = 1/2 plateau regions 
respectively. Orbital ordering may induce a further lowering 
of the symmetry |21j . 



effective Ising model 



(3) 



on a cubic lattice with nearest neighbor and next nearest 
neighbor antiferromagnetic interactions J\ = Vi/4 and 
J2 = ^/4, respectively, between Gi = ± 1 pseudospins. 
S = becomes a t = — 1 and S z = 1 becomes = + 1 via 
the transformation rii = (1 + c^)/2 and the adoption of 
an effective magnetic field ft = \q\ibH — \ A — |Vi — 3V2 
(shown schematically in Fig. [3]). We note the orbital 
orientation-dependent terms neglected in Eq.Q are rel- 
evant for deciding the orbital ordering that accompanies 
each spin ordering at low enough temperatures. 

Provided the nearest neighbor coupling dominates (i.e. 
Ji > 4J 2 ), the ft = ground state ft = is a two- 
sublattice Neel order, whose actual magnetization is half 
of the fully saturated value m = - corresponding to 
(<7i) = 0. We associate this ordering with the measured 
magnetization plateau in the region above ~ 75 T shown 



magnetization plateau m tne region above ~ 
in Figs. [5] and [l] Intermediate plateaux with 
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\ and 



m = I (in which (c^) = =p|) become viable for finite 
J2. On transforming back to the original Hamiltonian 
(considering both positive and negative ft), we therefore 
realize a cascade of first order transitions between dif- 
ferent phases with critical fields — H c \, H C 2, H c s and 
i^ C 4 in Fig. [3] — whose relative separations depend on the 
nearest and next nearest neighbor repulsions V\ and V2. 

Fig.|3]shows quantitative consistency between the mea- 
sured magnetization steps and those of the model owing 
to the quenching of the orbital contribution to the g- 
factor (such that g ~ 2) by a JT distortion of the 06 oc- 
tahedra encasing the S z = +1 sites. On associating the 
first two transitions (H c i and iJ C 2) with the two transi- 
tions observed experimentally, we estimate V2 ~ 1.2 K. 
The lack of further transitions in motor-generator-driven 
magnetic fields extending to 97 T (in Fig. [I]), and in a 
single turn magnet system delivering fields to ~ 140 T 
(shown in the inset to Fig. [TJb)), implies that V\ > 25 K. 
On enumerating V\ over the six nearest neighbor Co 
atoms in the cubic perovskite, we arrive at a lower limit 
< 150 K for the energy scale over which collective effects 
influence the phase diagram of LCO. 

Since the S z = +1 ions form a bipartite BCC struc- 
ture, the dominant antiferro-orbital interaction gives rise 
to a staggered orbital order in the m = 1/4 plateau. 
On the other hand, antiferro-orbital interactions are frus- 
trated in the second m = 1/2 plateau as the S z = +1 
ions form a non-bipartite FCC lattice. Nonetheless, our 
Monte Carlo simulations found a partial orbital ordering 
with a layered structure. 

The different crystalline structures in Fig. j3^b) (e.g. 
BCC and FCC) anticipated for each of the phases of 
LaCoOa provide an explanation for the discontinuous 
changes in the lattice and hysteretic behavior. Each mag- 
netization plateau in Fig. [3ja) corresponds to a different 
optimal sublattice arrangement of S = and S z = +1 
spin states. Orbitals in the S z = +1 sublattice are more 
voluminous and directional in nature, causing the lattice 
to expand in response to the pressure exerted by their 
increased density (i.e. V\ and V2). The reduced expan- 
sion at the second transition (^2) suggests the predom- 
inantly repulsive interaction is at least partially compen- 
sated by an attractive antiferro-orbital interaction, which 
is able to afford a more efficient arrangement of the or- 
bitals at higher densities [2T] , 

Finally, we reiterate that the spin value S of the first 
excited multiplet is still matter of debate [I5j |30l I3T] . 
While the generic form of the effective Ising Hamilto- 
nian of Eq. ([3| does not depend on the value of S, one 
would have to include further neighbor repulsions among 
magnetic configurations in order to explain a ~ 0.5/i^ 
magnetization plateau for S = 2 (~ 0.5/i# is 1/8 of the 
saturated value in this case). Therefore, our original as- 
sumption of S = 1 leads to a simpler and more natural 
explanation of the measured plateaux. We also note that 
the orbital physics and magnetostrictive properties are 
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likely to be very different for S = 2 owing to the e g or- 
bit als no longer being JT active. 

The experimental evidence therefore suggests that col- 
lective behavior involving two spin states leads to dif- 
ferent crystalline arrangements that can be tuned with 
magnetic field. While magnetostriction associated with 
field-tuned orbital order has been reported in mangan- 
ites [27], the present behavior in LaCo03 is different in 
that a strong coupling of spin, orbital and lattice degrees 
of freedom occurs in a Mott insulator. An entirely differ- 
ent type of functionality results: multiple field-tuned (dif- 
fusionless) transitions giving rise to a rapidly switchable 
strain. Spin state crystallization is likely to be a general 
property of crystalline materials subject to a spin state 
transition under extreme conditions [3] [4[ 8, 9 , causing 
such materials to potentially become metastable in the 
vicinity of a phase transition and vulnerable to a sudden 
release of mechanical energy. 
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SUPPLEMENTARY INFORMATION 
Further Experimental details 

Magnetization 

Magnetization to « 63 T is measured in a 'short-pulse' 
capacitor-driven pulsed magnet over temperatures rang- 
ing between 1.7 and 55 K using an extraction magne- 
tometer [20] consisting of concentrically arranged coun- 
terwound coils of ~ 1000 turns. Magnetization studies to 
~ 97 T are performed in a generator-driven nondestruc- 
tive magnet at Los Alamos National Laboratory [19]. 
The magnetic susceptibility is measured inductively as 
an induced voltage, with the discontinuous nature of the 
transitions typically yielding a signal in the form of a 
delta function of height ~ 10 V at each of the transi- 
tions. Susceptometer measurements in the 140 T single- 
turn magnet system [32] were made using a compensated 
coil consisting of a few turns. 

Magnetostriction 

Magnetostriction measurements are made using a 
fiber-optic Bragg grating (FBG) attached to the sam- 
ple using cyanoacrylate bond [33] . The single-mode fiber 
contains a 1 mm long Bragg grating and is manufactured 
by SmartFibre, UK. The reflection of light at the Bragg 



wavelength shifts when the sample length changes. Op- 
tical spectra are collected with a fast 1024 pixel InGaAs 
infrared line array camera, which is read out at 47 kHz. 
Mounted on a fiber-coupled research spectrometer of fo- 
cal length 500 mm with a 600 groove/mm -1 diffraction 
grating, the 25 jam pixel pitch results in a dispersion, 
SX of 0.056 nm per pixel at 1550 nm. To illuminate the 
FBG, we use a super luminescent diode broadband source 
with an output power of 30 mW in the wavelength range 
1525-1565 nm. 

Dissipation effects 

The hysteresis at each of the transitions evidences 
dissipation resulting in the unavoidable generation of 
heat during the experiments. This irreversible heating is 
likely partially compensated by cooling due to the mag- 
netocaloric effect, which in turn is caused by the field- 
induced closing of the spin gap. The sample is partially 
able to impart the resultant heat to the surrounding cryo- 
gen medium in times of order 1 ms at base temperature 
(T « 1.6 K) where the sample is immersed in superfluid 
4 He. Owing to the rapid succession of the two transitions 
on the rising magnetic field where the rate of change of 
magnetic field is highest, not all of the heat generated 
by the first transition is able to escape the sample be- 
fore the second transition, possibly contributing to the 
temperature-rounding of the second transition on the ris- 
ing magnetic field. The slower ramp rate of the falling 
magnetic field enables more heat to dissipate to the cryo- 
gen between the transitions, enabling both transitions to 
be clearly observed with reduced rounding. 

Maxwell relation and field-induced volume expansion 

The Maxwell relation dV/dH\ T , P = -dM/dP\ T , H , 
where V is the sample volume, H is the applied magnetic 
field, M is the sample magnetization, and P is the ap- 
plied hydrostatic pressure, provides a means to check our 
data for self consistency. We know from prior results [34] 
that dM/dP < for temperatures T < 100if , caused by 
the increase of the crystal field gap with pressure. Hence, 
an increase in the sample volume with magnetic field is 
expected. This is consistent with our experimental obser- 
vation of the longitudinal magnetostriction in high fields, 
and what we expect from our model. 

Details of the Model 

Background 

In the model, we assume the single-ion ground state 
of Co 3+ ion to be non-magnetic (S = 0) while the first 
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FIG. 4: Supplementary Information: A schematic rep- 
resentation of higher order interaction terms. Anti- 
ferro and ferro-orbital configurations on a nearest-neighbor 
bond along z-axis (a and b) and a next-nearest-neighbor bond 
in the xy plane (c and d). The corresponding interaction en- 
ergies are Vi, V/ , V?, and V/, respectively, where the sub- 
script '1' refers to nearest neighbor and '2' refers to second 
nearest neighbor. 



excited state is magnetic with 5=1. The latter has 
an orbital degree of freedom since there is one electron 
shared between two e g orbit als. We further assume that 
the energy gap (A) between the 5 = 1 and nonmagnetic 
configurations is large in comparison to the exchange in- 
teractions between different magnetic configurations of 
Co 3+ ions. Applying a magnetic field H lowers the energy 
of the doubly degenerate S z = 1 triplet states (the e g or- 
bitals being d% z 2_ r 2 and d x i_ y i) such that their energy 
becomes comparable to the 5 = nonmagnetic state. 
A low-energy effective model is obtained by projecting 
the original Hamiltonian into the subspace of the three 
low-energy states. 



Relevant terms in the Hamiltonian 

A single 5 = 1 ion embedded in the host matrix ac- 
quires additional energy via JT coupling with the lo- 
cal octahedron and exchange interaction with its non- 
magnetic neighbors, amounting to a slight modification 
(A) of the energy gap (A) entering into the single ion 
component (Hi) of the Hamiltonian given by Eq. (1) in 
the main text. As the density (rii) of 5 = 1 states in- 
creases, the interactions between them leads to the for- 
mation of spin-state superstructures and orbital order, 
necessitating a second component of the Hamiltonian 



i,3 



mechanism is the elastic interaction mediated by long- 
wavelength phonons. Short-range interactions are usu- 
ally orbital dependent and have contributions from both 
superexchange and JT couplings. 

The origin of the isotropic (Vij) term is discussed in the 
main text. The orbital-dependent interaction potential in 
Eq. Q consists of two terms: 



ViAsi,Sj) 



V a - o a 

IJ Si ,Sj 



(5) 



where the ferro- and antiferro-orbital bond operators are 
given by 



Vij + VlAs^Sj) 



(4) 



°Ls 3 = (6) 

°l,s 3 = S aitaij (1 - 5 SjiOCi .) + (1 - 6 8itaij ) 5 Sj , aij . 

Since the orbital interaction is generally anisotropic, the 
variable ol^ = x,y,z indicates the orbital flavor (orien- 
tation) which has distinct symmetry properties with re- 
spect to the other two on a given bond (ij). For exam- 
ple, OLij = z (ds z 2- r 2 orbital) on a nearest-neighbor bond 
along z-axis or a next-nearest-neighbor bond in the xy 
plane. Orbital configurations in which the a-orbital is 
empty on both sides are neglected in Eq. ^ since such 
terms are irrelevant for orbital ordering. 

The orbital-dependent interaction potentials (V/j and 
Vf!-) result mostly from the anisotropic superexchange 
and JT interactions, as described in more detail in 
the following section. Since the elongation of 06 oc- 
tahedra is the source of local anisotropic stress, the 
resultant phonon-mediated interaction is also orbital- 
dependent [35] thereby contributing to both V/j and V^. 
The interaction potentials for two S z = +1 ions on the 
nearest neighbor and next-nearest neighbor bonds are 
summarized in Fig. [4j The electronic contribution to 
the interaction potentials becomes negligible beyond a 
few nearest neighbors as higher-order hoppings are re- 
quired for such virtual exchange processes. The orbital- 
dependent interaction potentials (V/j and V^) result 
mostly from the anisotropic superexchange and JT inter- 
actions, as described in more detail in the following sec- 
tion. Since the elongation of 06 octahedra is the source 
of local anisotropic stress, the resultant phonon-mediated 
interaction is also orbital-dependent [35] thereby con- 
tributing to both V/j and Vgj. The interaction poten- 
tials for two S z = +1 ions on the nearest neighbor and 
next-nearest neighbor bonds are summarized in Fig. [4] 
The electronic contribution to the interaction potentials 
becomes negligible beyond a few nearest neighbors as 
higher-order hoppings are required for such virtual ex- 
change processes. 



The effective Hamiltonian 



describing interactions between different sites. H2 is 
the same as Eq. (2) in the main text. The dominant 



While nearest-neighbor ferro- and antiferro-orbital in- 
teractions are expected to be significant, these terms do 
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not play a role in the formation of the first two field- 
induced phases. The density of S z = +1 orbital states 
is too low to support the existence of such bonds. It is 
therefore reasonable to assume that the net interaction 
is repulsive, with larger orbitals tending to repel each 
other, enabling us to construct an effective Ising model 
Hamiltonian (% e ff) consisting of nearest neighbor (Vi) 
and next-nearest neighbor (V2) repulsive terms as given 
by Eq. (3) in the main text. 

The magnetization curve m(H) of the original model 
is shown in Fig. 3a in the main text. The width of the 
m = 1/4 and 3/4 plateaus is AH1/4 = AH 3 ^ — 24 J2 = 
6V2, indicating that both states are stabilized by second- 
nearest neighbor interactions. 



find 



V 2 a 



SjTt 2 



dda 



Udd(Udd -<5jt)' 



(8) 



where tdda is the direct dda hopping between t 2g orbitals 
and Udd is the on-site Coulomb repulsion of d electrons. 
Unlike the case of m = 1/4 plateau, The antiferro-orbital 
interaction V 2 a remains frustrated as the S z = 1 ions in 
the half-plateau form a non-bipartite face-centered cubic 
(FCC) lattice. Nonetheless, our Monte Carlo simulations 
found a partial orbital ordering with a layered structure. 
Detailed studies of the orbital order will be presented 
elsewhere. 



Orbital ordering 

In the m = 1/4 plateau phase, likely corresponding to 
the first field-induced phase (SSC1) in Fig. 3b in the main 
text, the closest distance between S z = +1 ions is r = 
\^3a (i.e. third- nearest neighbors). Since superexchange 
and JT interactions are negligible at this distance, we 
resort to the anisotropic elastic interaction discussed in 
Ref. 35 to calculate the effective interaction potentials: 



V f (r) = 3cr r/87rr 3 , 
V a (r) = -3cr r/167rr 3 , 



(7) 



where T = (3cn + 4c^)/{cn + 2C44), en and C44 are 
elastic moduli of the crystal, and <Jo is the strength of 
the local stress. Since the third-neighbor bonds are par- 
allel to the (111) directions, the three 3l 2 — r 2 orbitals 
(I = x, y, z) are equivalent with respect to the bond. The 
corresponding bond operators are simply 0{. s . = S SiySj 
and Og.^ s . = 1 — S SijSj . The elastic interaction given by 
Equation ^ shows an additional energy cost for ferro- 
orbital configurations, we therefore expect the S z = +1 
ions at the m = 1/4 plateau to form a bipartite body- 
centered cubic (BCC) lattice with a doubled lattice con- 
stant and a staggered orbital order. Since there are three 
possible orbital orientations, there are six different stag- 
gered orbital configurations that have the same energy 
in absence of spin-orbit coupling: d 3zy 2_ r 2, <i 3zy /2_ r2 with 
v, v' = {x, y, z} and v ^ v' . 

In the m = 1/2 plateau phase associated with the sec- 
ond field-induced phase (SSC2) in Fig. 3b in the main 
text, the next-nearest neighbor antiferro-orbital compo- 
nent (V<2 in Fig.[4| is expected to be relevant, which has 
a contribution from superexchange. While indirect hop- 
ping via O ions only contributes to sixth order to the ex- 
change energy between next-nearest-neighbor Co 3+ ions, 
direct hopping between t 2g orbitals contributes to second- 
order. On taking into account the splitting Jjt of t 2g 
levels due to the elongation of the local octahedron, we 



Calculation of bond energy 

Exchange interaction 

Here we compute the exchange energies for a pair of 
Co 3+ ions on a nearest neighbor bond. Figure [5] shows 
six inequivalent configurations of Co ions on a nearest- 
neighbor bond along the z-axis. The directional depen- 
dence of e g orbitals results in a highly anisotropic electron 
hopping, which in turn leads to an anisotropic exchange 
interaction. To simplify the calculation, we consider only 
the dominant pda hopping between e g and p orbitals of 
the 2 ~ ions. The resultant exchange energy of a z-bond 
thus depends on whether the d 3z 2_ r 2 orbitals are occu- 
pied. We expect E a = = Ed = cq< E c = E e = e\ and 
Ef = 62- Here the subscripts of e indicate the number 
of occupied d 3z 2_ r 2 orbitals on the bond. For simplicity, 
we also neglect Hund's coupling. 

We employ standard degenerate perturbation theory 
to calculate the bond energies up to fourth order in t p d- 
For configurations a, b and d in Fig. [5j the second and 
fourth order corrections to the bond energy are 



and 



S 2 ) 



4t 2 , 



(A p d — U pp ) ' 



(A, 



U ) 2 

Upp J 



(Apd — Up 



(9) 



(10) 



(2Ap d — Upp) (2A pd + Udd - Up 



Here, U pp and Udd are the Coulomb repulsion energies 
for double occupation of the p z and d 3z 2_ r 2 orbitals re- 
spectively, while A p d is the diagonal energy difference 
between the two orbitals. We have neglected corrections 
of order Jh/U where Jh is the Hund's exchange cou- 
pling. These perturbative corrections to the energy can 
be understood to result from the virtual hopping of the 



FIG. 5: Supplementary Information: A schematic showing the six inequivalent relative configurations for a 
pair of Co 3+ ions. Panel a shows a pair of nonmagnetic S = ions, b and c show a pair of S = and S = 1 ions with 
d x 2_ y 2 and d 3z 2_ r 2 orbitals, respectively, d, e, and f show a pair of S z = 1 ions with orbital configurations (d x 2_ y 2 , d x 2_ y 2), 
(d 3z 2_ r 2 , d x 2_ y 2 ), and (d 3z 2_ r 2 , d 3z 2_ r 2), respectively. 



p electrons via the ion to the neighboring Co sites and 
back. Similarly, we obtain the energy corrections 



ol pd 



(A pd — U pp ) ' 



(ii) 



where 7 is a coupling constant. By integrating out the 
displacement u, the energy gain for the six different con- 
figurations are 



E h 



11 
8fc' 



E, 



71 

2fc' 



E, 



97^ 
" 8fc ' 



(Apd - Upp) 2 (2A pd + Udd — Upp) 

Kd 

(A pd -Upp) 2 (2Ap d -Upp) 



L pd 



9t 4 



(A pd - Upp) 2 U dd (A pd - Upp) 3 



M) - 



it 2 



(A pd — Upp) ' 

(A pd - U pp ) 3 ■ 



(12) 



(13) 
(14) 



Jahn-Teller interaction 

We compute the energy gain of a bond via Jahn-Teller 
coupling to the 06 octahedra. For simplicity, we adopt 
the so-called Einstein phonon model for the displacement 
u of the ion, namely E e = k\u\ 2 /2. For a bond oriented 
along the z-axis, the JT couplings for the d 3 ^2_ r 2 and 



E a — E d — E 



f 



0. 



(15) 



Calculation of interaction potentials 



Here we derive the effective interaction potentials be- 
tween a pair of Co 3+ ions. In particular we will focus 
on the exchange contributions to the potential energy 
We introduce a pseudospin 1/2, r, to describe the or- 
bital degrees of freedom of the S z = +1 ions such that 
\r z = ±1/2) correspond to occupied orbital \3z 2 — r 2 ) 
and \x 2 — y 2 ), respectively As can be seen from the 
previous discussion, the exchange energy of a given bond 
depends only on the occupation of the active orbital (e.g. 
d% z 2_ r 2 on a z-bond). We thus introduce an operator P v 
for the three different nearest neighbor bonds (y = x, 
or z) which projects onto the relevant orbital state. On 
introducing the following unit vectors 



n 



x/y 



± — x z, 

2 2 ' 



the bond projection operators become P v — 
orbitals are ju z and —ju z /2, respectively [36], The exchange Hamiltonian is then given by 



(16) 



r • n, y . 
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v (ij)\\v 

+Gi [rii(l - raj) + (1 - ra^ra^] + G 2 [ra*(l - ra^) n • h u + (1 - ra^ra^ Tj • n^]|. (17) 



r 



Here i^i = e 2 + e - 2ei, K 2 = (e 2 - e )/2, K 3 = (e 2 + 
e + 2ei)/4, Gi = (ei + e )/2, and G 2 = (ci - e ). 

In the regime dominated by nonlinear JT interactions, 
quantum fluctuations of the orbitals are quenched by the 
lattice distortion. We thus replace the pseudospin oper- 
ator by a classical 0(2) vector of length |r| = 1/2. The 
local elongation of the 06 octahedron aligns the vector 
to one of the three directions h u . Each corresponds to 
the three Potts states discussed in the main text. Noting 
that t i • yi u = 1/2 for Si = v and —1/4 otherwise, the 
effective nearest neighbor interaction potentials in Equa- 
tion (El) can be easily read off from H ex as 



J_ 

16' 



3J 
16 ' 



f 



~16~' 



where the exchange constant J is given by 



J = e 2 + e - 2ei 



(18) 



(19) 



(A p d - Up 



A-U v 



U dd 



and the correction to the gap energy is 



A ex = 3(ei - e ) 



OL pd 



(A^ - Upp) 



+ 0{ti d ). (20) 
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